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Abstract 

The paper considers the motion of a body under the influence of gravity and drag of the 
surrounding fluid. Depending on the fluid mechanical regime, the drag force can exhibit a 
hnear, quadratic or even more general dependence on the velocity of the body relative to the 
fluid. The; case; of quadratic drag is substantially more complex than the linear case, as it 
nonlinearly couples both components of the momentum equation, and no explicit analytical 
solution is known for a general trajectory. After a detailed account of the literature, the 
paper provides such a solution in form of a series expansion. This result is discussed in 
detail and related to other approaches previously proposed. In particular, it is shown to 
yield certain approximate solutions proposed in the literature as limiting cases. The sohition 
technique employs a strategy to reduce systems of ordinary differential equations with a 
triangular dependence of the right-hand side on the vector of unknowns to a single equation 
in an auxiliary variable. For the particular case of quadratic drag, the auxiliary variable 
allows an interpretation in terms of canonical coordinates of motion within the framework 
of Hamiltonian mechanics. The proposed reduction strategy also permits to devise analytic 
solutions for more general drag laws, such as general power laws or power laws with an 
additional linear contribution. Furthermore, a generalization to variable velocity of the 
surrounding fluid is addressed by considering a linear velocity proflle, for which a solution 
in form of a series is provided as well. Throughout, the results obtained are illustrated by 
numerical examples. 



Corresponding author. Email address: Shouryya.Ray@mailbox.tu-dresden.de 



Contents 



1 Introduction [3] 

2 A reduction principle for certain coupled systems of ordinary differential 
equations [5] 

3 Projectile motion with quadratic resistance [7] 

3.1 Mathematical formulation [7] 

3.2 Solution E 

3.2.1 Reduction to a Single Scalar Equation |5] 

3.2.2 Solution of the Reduced Problem [TU] 

3.2.3 Relation to the Hodograph Solution [T2] 

3.3 Properties [12] 

3.3.1 Integral of Motion [12] 

3.3.2 Limiting Behaviour [TJ] 

3.4 Numerical Illustration [E] 

4 Generalizations 1161 

4.1 Generalised Drag Laws [TU] 

4.1.1 Mathematical Problem [18] 

4.1.2 Solution M 

4.2 Accounting for a Velocity Profile [20] 

4.2.1 Equations of Motion [2D] 

4.2.2 Solution M 

4.2.3 Properties [52] 

5 Conclusions and Outlook 1241 
References 1271 



2 



1 Introduction 



Projectile motion constitutes a very elementary problem of classical mechanics; as such, occupying 
a central place in the works of Niccolo Tartaglia, Galileo Galilei and Sir Isaac Newton, to name 



but a few. Although Hardy (19401 scathingly remarked in his A Mathematician's Apology that 
the science of ballistics were "repulsively ugly and intolerably dull", he still admitted that it 
does demand "a quite elaborate technique". The latter statement perhaps explains why it was of 
interest to mathematicians such as Johann Bernoulli, Leonhard Euler, Adrien- Marie Le Gendre 
and Johann Heinrich Lambert, throughout the course of its long history. 

The basic physical problem consists of calculating the motion of an object moving in a 
fluid medium, experiencing the (constant) gravitational force of the Earth and the resistance 
of the fluid, where the term "fluid" encompasses both liquids and gases. In Book II of his 
Philosophi<£ Naturalis Principia Mathematica (1687), Newton gives the first rigorous treatment 
of this problem using mathematical physics. A simple force balance can be used to account for 
the two principal forces, namely gravity and drag, and, using his laws of motion (established in 
Book I of the same work), the equations of motion can be formulated in vectorial notation to 
read: 

mp^=Fa + F^. (1) 

Here, rup denotes the mass of the moving object, v the velocity vector of the object in a suitable 
inertial frame of reference with the initial condition v(t = 0) = 1^0 and Fq = rripQ is the gravity 
force, where g denotes the constant gravity acceleration vector. A Cartesian coordinate system 
(x, y) is usually chosen, with x the horizontal coordinate and y the vertical coordinate, such that 
g may be supposed to be oriented in the negative y direction, i.e. g = (0, —g). The equations 
of motion then constitute a system of ordinary differential equations in the components of 
V = (v^,Vy). Fp, denotes the drag force. The nature of the drag force then determines the nature 
of the equations of motion and thereby the mathematical problem at hand. The statement of 
the solution is completed by imposing the initial conditions v(t = 0) = Vg. 



Newton (1687) considered three forms of drag, namely linear, quadratic and a superposition 



of the two. If the fluid resistance varies in direct proportion with the relative velocity (i.e. linear 
dra^ , the differential equations for and Vy decouple due to the linearity of the drag force and 
constitute a system of linear first-order ordinary differential equations. The solution can be readily 
found using elementary methods (cf. textbooks on the theory of ordinary differential equations. 



e.g. Walter [2000 §2.1). Should the resistance, however, vary as the square of the velocity 
(quadratic dra^, one has to deal with a system of nonlinearly coupled ordinary differential 
equations, and the solution necessitates a more involved approach. Newton himself was unable 
to solve the problem, but his contemporary, Johann Bernoulli, solved it after being challenged 



^This is called Stokes drag after Stokes (1851), who derived the force exerted on a moving particle in a viscous 
fluid for vanishing Reynolds numbers. 

Unlike linear drag, for which Newton was unable to give any physical explanation (this was addressed much 
later by Sir George Gabriel Stokes in 1851), the quadratic drag was explained heuristically in the Principia in a 
manner much closer to the modern textbook approach. Such a law is valid for Reynolds numbers roughly between 
10"^ and 10^, for example. 
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by the British astronomer John Keih (Bernouhi 1719 1 . BernouUi's solution parameterizes the 



absolute value of the velocity over the trajectory slope (angle) — and is, hence, implicit. This 
solution is also known as the hodograph solution. The hodograph, being the locus of the tip of the 
velocity vector with the other end held fixed, has precisely the two afore-mentioned quantities 
as polar coordinates. Furthermore, the solution contains quadratures which must be evaluated 
numerically. Despite these drawbacks, it is the standard and by far most widely cited solution. 



cf. Synge and Griffith (19491, Hayen (20031, Benacka (20101 and Benacka (2011 1 . 



Much of the literature on projectile motion with quadratic drag after Bernoulli comprises 
efforts to calculate approximate solutions based on the hodograph solution. Various exact and 
approximate implicit formul^^ using miscellaneous series and approximation techniques can be 
found in the extensive literature available on this subject, the most well-known of which are due 



to Euler (1745), Lambert (1767), de Borda (1772) and Le Gendre (1782). The exposition of their 



complete results is beyond the scope of this essay. The interested reader is referred to Isidore 
Didion's Traite de balistique, sect. V, § III-IV, where these results have been painstakingly 



compiled. A Review on the historical aspect of the problem is also provided by Hackborn (2006). 



From the late 19th century onwards, numerical and empirical methods superseded the more 
analytical approach of Bernoulli and Euler. Nevertheless, the problem continues to afford scope 
for an instructive application of mathematical analysis to physics. An approximate formula 



of the trajectory for low quasi- horizontal paths is found in Lamb's Dynamics (1923). [Parker 



(1977) rederived the result and provided — to the best of our knowledge — for the first time. 



approximate explicit solutions for fiat trajectories and short duration. He also arrived at the 
results of Bernoulli by following a considerably different set of variable transformations. Another 



recent approximation is due to Tsuboi (1996), who considered perturbative solutions of the first 



order for small angles of release. Lastly, an explicit exact albeit semi- analytical solution was 



given by Yabushita et al. ( 2007 ) using the recently developed homotopy analysis method of Liao 



(2004). 



Summa summarum, it would appear that, in spite of the long history and the substantial 
amount of material available on this subject, an explicit exact analytical solution cannot, to the 
best of our knowledge, be found in the existing literature. Providing such a solution would hence 
be of interest. 

Drag laws more general than linear or quadratic have, for obvious reasons, received significantly 



less attention. Bernoulli (1719) provides an exact implicit solution using the hodograph technique 



for cases where the drag force varies as a general power of the velocity. Such laws are valid for 
high velocity subsonic and some cases of supersonic projectile motion, with Reynolds numbers 



substantially higher than 10 (Weinacht et al. 2005). For Newton's generalisation involving a 



superposition of quadratic and linear drag, no exact solutions can, to the best of our knowledge. 



be found in the literature. According to Clift et al. (1978), it may be instructive to consider even 



The expressions derived therein are imphcit in the sense that the unknown quantities (the velocity components 
and Vy or alternatively the cooridinates of the position of the particle x and y) , instead of being expressed 
as functions of time, are parameterized over some other auxiliary quantities, for example trajectory slope or 



slope-angle. Of particular interest is Lambert (17671, where y and t are parameterized over x 
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more general forms, especially for motion with Reynolds number between 1 and 10^. Although of 
no significance to exterior ballistics, projectile motion in the transition regime between Newtonian 
and Stokesian drag may be studied in order to gain general understanding of particle motion in 



sediment transport phenomena — somewhat along the lines of Nalpanis et al. (1993). 

A final generalisation that we shall consider in this paper deals with the effects of a velocity 
profile on projectile motion, as no studies investigating the influence of a velocity profile of any 
kind on projectile motion with nonlinear drag could be found in the literature. Although it 
would be necessary, for the analysis to be exhaustive, to consider general drag laws and profiles, 
we restrict ourself to the quadratic drag law, for it is, in some sense, the simplest nonlinear drag 
law, and linear velocity profiles. 

The paper is structured as follows: We begin by establishing a useful principle for the 
reduction of a certain class of nonlinear systems of ordinary differential equations to a single 
equation. It is then exploited to solve the basic problem of projectile motion with quadratic 
resistance. The various inter-relationships between the explicit analytical solution thus obtained 
and the previous attempts are then examined. Subsequently, the limiting behaviour of the 
solution in cases of physical interest is investigated. An integral of motion is derived, its relation 
to historical forms found in the literature is elucidated and a physical interpretation is suggested 
based on the present form. 

In the next part, two possible generalisations of the basic problem are proposed. First, the 
approach employed in solving the quadratic drag problem is used to tackle more general drag 
laws. Thereafter, the scenario is generalised by including the effect of velocity profile and the 
governing equations are solved by virtue of the reduction principle established earlier. Finally, 
some interesting properties of the solution are discussed. 



2 A reduction principle for certain coupled systems of ordinary 
differential equations 

Let X = (X-^, X2, X^, . . .), be an n-dimensional vector- valued function of t G R and let it be 
determined by the initial value problem 

dX 



dt 



+ PX = Q, X{t = 0) = Xo 



(2) 



Here, P : II — )■ R denotes a scalar real- valued function and Q = (Qi, Q21 Qi, . . .) : H — )■ E." is 
a vector-valued function. Were P and Q known functions of t, the equations would constitute 
a linear system of ordinary differential equations and there would be no difficulty in applying 
the method of variation of parameters (Lagrange 1809a|b , 1810) in order to write down the 
solution for each and every component of Eq. ^ in terms of P, Q and appropriate quadratures 
involving the two, thus 

Xo + ^ Qexp(^^ Pds^ dr 



X 



exp 



PdT 



(3) 
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with Xq denoting the vector of initial values. 

For the purpose of this investigation, however, it is essential, that we consider a more 
general class of equations. Assume now, that P is a scalar function of X and (possibly) t, i.e. 
P = P{X] t). Furthermore, assume that Q = Q{X; t), with the constraint that Vj^Q, which is a 
tensor of rank two, be a strictly lower triangular matrix, i.e. 

Qi = Qi{t) 
Q2 = Q'zi^ij t) 

Q, = Q,{X,,X^-t) (4) 



In that case, variation of parameters is inadequate and the expression in Eq. ^ can hardly be 
called a solution. Ignoring for a moment that the equations are no longer linear and naively 
using variation of parameters on the first component of Eq. ^ , the following formal expression 
is obtained: 

^1.0+ fQ,{T)exp(rP{X;s)ds) dr 

X, = ^ ^ . (5) 

exp(^y^ P(X;r)drj 

The only unknown quantity is the exponential of / Pdt, and it will later be clear that it is a key 
quantity in this method. Therefore, let 



exp(^^V(X;r)dT) (6) 



for the sake of brevity. Now, let us take the liberty of expressing Eq. ^ more succinctly by 
writing X-^ = fi{4>). This can be taken over to the next component of X yielding 

fQ,{X,;T)exp(rP{X;s)ds)dT X,,, + f Q.UMY.r) <pdT 
_ ■Iq \Jo / _ Jo ^^-j 



exp^V(X;T)dr^ 



This, again, is denoted as X2 = A ((/>)• Likewise, the i-th component of X then may be expressed 
recursively in a similar fashion. From a general point of view, this procedure defines an n-tuple 
of mappings / = (/i, /a, /g, ...):/—)■ D from the space I of positive integrable functions to the 
space D of differentiable functions, so that ry E / is mapped onto 



^1,0+ / Qi{r)vdT 
J a 



m ^ — '-^ ^ = 1 (8) 



^..0+ / Q,[/i(??),/2(r/),...,/.-i(r/);T]77dr 

fM = 1< i ^ n (9) 
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The solution of Eq. @ can now be formulated as 



X = /((/)) 



(10) 



with (j) defined according to Eq. Q . Inserting this into the first component of Eq. ^ gives 



(11) 



dt 



In this equation, the only unknown quantity involved is (p. Once (j) is known, X is readily 



system of ordinary differential equations has been reduced to a single equation in an auxiliary 
variable appropriately defined for that purpose. For ease of reference, we call it the resolvent 
variable and the equation governing the auxiliary quantity the resolvent equation. 

Remark on scope of application. The above formalism relieves one from the trouble of solving 
a whole set of coupled ordinary differential equations and rather allows one to concentrate the 
investigation on a single equation for a scalar quantity only. Since / involves quadratures over 
(possibly) non-trivial kernels, the resolvent equation is, in general, an integro-differential equation 
and, therefore, may not be trivial itself. It clearly depends to a large extent on the nature of P 
and Q as to how difficult the new equation is. It will, however, be seen in the following that in 
the case of the particular class of problems considered here, the formalism may be exploited 
to obtain elegant analytical solutions of the respective original systems of ordinary differential 
equations. 

3 Projectile motion with quadratic resistance 
3.1 Mathematical formulation 

The general force balance ([T]) discussed above is now supplemented by specifying the drag force. 
The quadratic resistance law can be written in the form 



Here, Q, is the constant drag coefficient characteristic of the projectile's shape, A is the cross- 
sectional area and g is the density of the fluid. The velocity of the fluid which constitutes the 
resistive medium is denoted by V, which, at present, we assume to be constant in space and 



calculated using Eqs. (|8])-(10). Thus, the original problem involving a non-linearly coupled 





time. Equation ([T]), therefore, reads 



dv 
lit 



a \\v — V\\ {v — V) + g , 



(13) 



where the constant a is a shorthand for ^ 



C^gA/mp. Finally, we introduce the quantity 



u = V — V , 



(14) 
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which is the relative velocity between the particle and the fluid. Since V is a constant vector, 
this yields 

— = —a\\u\\u + g. (15) 
or, when using the components of u = {u^,Uy) 



du 



" +aUyJul +ul 



-9- 



(16) 
(17) 



Without loss of generality, the coordinate system may be placed at the starting point of trajectory. 
The initial conditions are given by 'u(O) = t'(O) — V , or, in components. 



(18) 
(19) 



Observe that a purely vertical motion (in the reference frame moving at the fluid velocity) is 
obtained if u^ q = 0. That case, again, is described by one component only and the ordinary 
differential equation can be solved using separation of variables. Here, however, we are interested 
in the more general case and eliminate this situation by requiring that 



^^.,0 > 0. 



(20) 



With this condition, Eqs. (16|-(19| constitute the initial value problem to be solved in the 
following. 



3.2 Solution 

3.2.1 Reduction to a Single Scalar Equation 

Review of Parker's Approach. Since the present method of solution is somewhat akin in 



spirit to the one employed by Parker (19771, it is instructive to study the latter now. Division 



of Eq. (17) by (16), which is well-deflned due to the condition (20), yields, after appropriate 
algebraic manipulation, 

dUy/dt + g ^ Uy 
du^/dt 

This can be simplified using the quotient rule of differential calculus and subsequent separation 
of variables to give 

^ ^' (22) 



d( 



-g 



Upon integration and rearranging of terms, one arrives at the relation 

dr\ 



9 



(23) 
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Inserting this into Eq. ( 16 ) and substituting 



w = gf^-'^ (24) 
Jo u^fi 

yields the equation 

-^=ag^J\^w'. (25) 

This is the resolvent equation by Parker's method, i.e. the task of solving the original system is 
now reduced to solving this scalar equation subject to the initial conditions w[^) = —Uy ^/u^ Q 
and dtt;(0)/dt = gju^ ^. 



Modified approach. One may recast Eqs. (16) and (17) into the form considered in Section [2] 
by writing 

' Q = (0, -g) (26) 



Since X^Q = O2.2 is a null matrix, the condition of applicability is trivially fulfilled. It is, therefore, 
natural to chose 

(27) 



exp yuf+n^ dr 



in order to express u as 



u 



x,0 



\{u,,-gl^dr). 



(28) 
(29) 



Inserting into Eq. (16) yields 



dt 
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>dr 



(30) 



which is, as expected, an integro-differential equation. However, since the only kernel involved is 
unity itself, the equation can be reduced to an ordinary differential equation by rewriting it in 
terms of the quantity 



^ = 



'dr 



to yield 



d^ 



(31) 



(32) 



The initial conditions follow directly from the definition of the quantities ^ and (p as ^^(0) = 
and d^(0)/dt = 0(0) = 1. This is the new initial value problem equivalent to the original system 
which will be investigated further in the subsequent part of the paper. 

In view of the derivations made so far, two comments are now in order: 
(i) The result of Parker's approach and the present one are essentially equivalent, which is 
readily shown by observing that w = {g'P — Uyg) /u^q. Neither the form involving $ nor the one 
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involving w has a distinct advantage over the other. While the right hand side of the equation 



for w, (25), is somewhat simpler than that of Eq. (32), the initial conditions of the latter are 
slightly simpler. 

{ii) The difference between Parker's approach and the present one resides rather in the method 
than in the result. The present method can be readily generalised to tackle more involved 
problems, such as the one involving velocity profiles, as shall be shown later, while Parker's 
method is not so flexible. To be precise, Parker's method is only applicable as long as V„Q is 
a null matrix, whereas for the present method, it only needs to be a strictly lower triangular 
matrix. Therefore, this formalism may also be regarded as a means of embedding the solution 
techniques employed for the present problem into a more general framework. 



3.2.2 Solution of the Reduced Problem 



Since the right hand side of Eq. (32 ) can be developed in a power series around the origin {t = 0), 
the classical theory of ordinary differential equations then guarantees the existence of a power 
series expansion for ^ around t = 0, i.e. ^ = Og + a^t + Oa^^ + a^t^ + • • • . with 



1 d"^(0) 



(33) 



" nl dt" 

according to Taylor's theorem. The first two coefficients are determined using the initial conditions 



yielding Qq = Oanda^ = 1. Evaluating Eq. (32) at t = provides a direct relation for the third 
coefficient: 



a 



(34) 



Differentiating Eq. ( 32 ) once at t = and algebraic rearrangement gives 

aa = - (35) 

This process, continued ad infinitum, yields all the coeffients of the series. In general, differenti- 



ating both sides of Eq. ( 32 ) n times with respect to t gives 



a 



(n + 2)! dr 



(36) 



This, in conjunction with the initial conditions, constitutes a recursive formula for the coefficients 
of the power series expansion of ^, which solves the problem. For the sake of convenience, 
let Mg = ||wo|| = \ uifi + u^y be the Euclidean norm of the initial velocity vector and let 



arctan (Uj^ q/m^_q) be the initial angle of the trajectory. Using — u^ q /{d<P/dt) and 



10 



Uy = {Uy Q — g'P) /{d'P/dt), the final solution is as follows: 



^ = t + — - e 
1-2 



1 • 2 • 3 



gff^^o ^ cos^6'o - a^gup sin 
1-2-3-4 

Sagf^UQ ^ sin 6*0 cos^^q + a^g^ (sin cos^^q + sin^6'o + 3 cos''6'o) 
^ 1-2-3-4-5 



(37) 



05 sin 6*0 2 , aff^no'cos^6'o - a'5r?ioSin6lo 
1 + oMq t : — - — t H : — - — t + ■ ■ ■ 



1 • 2 



1 • 2 • 3 



u. 



05 sin 6*0 2 , ag^Uo^ cos' e^- a^guosm 00 

1 + OMq I ; t H t + ■ ■ ■ 



1 • 2 



1 + 



gt X 



OUq 



05 Sm fg 

1-2-3 



1 • 2 • 3 

ag'ug^ cos'Og — a^gu^ sin 3 



1 • 2 • 3 -4 



ag sin 6*0 oc/'^Mo ' cos^6'o - o'c/tXo sin „ 
1 + oMq i ; — — t H ; — — j;^ t + 



1 • 2 



1 • 2 • 3 



(38) 



(39) 



Finally, one can return to the corresponding absolute quantities with = + and v 



Uy + Vy. 



Remark on structure of solution. It is apparent that, upon truncating the power series of 
^ at n terms, one has a ratio of polynomials, where the numerator is of n-th order and the 
denominator of order (n — 1). The lowest neglected order in t is, therefore, n, so that it would 
be appropriate to refer to such an expression as an approximate of order (n — 1). 



Remark on Parker's implicit solution. Parker (1977) integrated Eq. (25) once to obtain 



dw 
df 



agJ C — ag (^w\/ 1 + w'^ + arsinh vuj , 



(40) 



where C is a constant completely determined by the initial conditions. Separation of variables 
yields 

r (41) 



C — ag + u>' + arsinh 



Perhaps in part due to the daunting nature of the integrand, he remarked that "even if the 
indefinite integral could be evaluated, we would not be able to invert the resulting expression to 
obtain an explicit solution". While it is probably inevitable that the integral cannot be evaluated 



analytically, it is interesting to note that due to a recent result from Dominici ( 2003 ) , it is now 



possible to obtain an explicit solution in the form of a power series from Eq. (41 ) without having 



to evaluate the integral in the first place. The main result of the cited reference states that if a 
function Y is implicitly given as 

^ dy 

X = i (42) 
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the relation can be solved for Y and has the form 



Y = Y{X = 0) + f{Y= 0) ^ = 0) 



X 



n + l 



(43) 



(n + l)!, 

where 2)"[/](y) denotes the n-th nested derivative of / with respect to Y, defined recursively as 

(44) 
(45) 



D°[/](y) = 

d{/.D"[/](y)} 



dY 



Writing, for the sake of convenience, w = dw/dt and applying Dominici's theorem yields 



w 



Uy,0 , 9 



^ + -^J2^[w]iw = -n,.oK,o)— — . (46) 
Due to the fact that w = {g^ — u^ ^) /u^ ^, the following alternative expression for ^ is derived: 



n! 



(47) 



n + l. 

Calculating the individual terms of this series and comparing them with the ones predicted by 



Eq. (37) shows that the two expressions are identical. 



3.2.3 Relation to the Hodograph Solution 

The hodograph solution due to Johann BernoullQ reduces the original system to a single ordinary 
differential equation with u = \\u\\ = Jul + Uy, 



yul + Uy, the Euclidean norm of the relative velocity 
vector, as dependent variable and 9 defined by tan 6 = dy'/dx' as the independent variable. 



Bernoulli's approach can be recovered from the present one as follows. Dividing Eq. (29) by (28), 
rearranging and using parametric differentiation yields the identity 



"^ = ^-°-^'^ =tang. 



(48) 



Together with Eq. ( 28 ) , this yields 



\/ulo + {uyo -g<py 



(49) 



This relation, along with the identity 



d^^ 




rd^\ 


d^ 




^d<P\ 


w ~ 


-( 

dt^ 


vd^J 


~ d^ 


-( 

d^^ 


. d^j 



d^ 



(50) 



It is difficult to know the exact method used by BernouUi because, as was usual during his time, he only 
provided the final result without any intermediate steps. All knowledge of his method, therefore, is based on the 



writings of his student Euler ( 1745 I, who appears to cite Bernoulli's results and extended them. 
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and an elementary (but cumbersome) application of the chain rule of differential calculus yields 



dn ^ au^ , ^, 

— =utan9 + -, 51 

dt^ g cos 8 

a so-called Bernoullj^ differential equation for u in terms of 9. Although the auxiliary equation 
here has an analytical solution in terms of elementary functions, in order to change from the 
parameterisation by 9 to the one by t, quadratures which cannot be evaluated analytically are 
necessary. This limits the usefulness of the approach. Nevertheless, the fact that one approach 
can be recovered from the other may be regarded as a consistency check. 



3.3 Properties 



Here we shall illustrate what information may be extracted from the analytical solution (37)-(39), 
with the ultimate goal of obtaining further insight into the properties of the motion. 



3.3.1 Integral of Motion 

Contrary to physical intuition, projectile motion under quadratic drag obeys a conservation law. 



This may be seen when using Eq. (50) to rewrite Eq. (32) as 



<^d(/--ad^V^<o + Ko-5^) =0, (52) 
which is readily integrated to yield 

^(/>' + C/.(^) = const. (53) 
Here, U,{<P) = —a J \J uI q + {Uy g — g^Y d^. The quadrature involved is elementary. An earlier 



and more common form for such an integral of motion is (Didion 1860) 

1 1 



2 a a + ^(^) = const. (54) 



2 cos' 9 sin 9 



where m ^ J\se^9\ d9. While the two forms are mathematically equivalent, Eq. @ may 
have certain advantages insofar as physical interpretation is concerned. Indeed, multiplying Eq. 



(53) with the mass of the particle, nip and defining nip U, = U, one obtains 

^TTip + U{<P) = const. (55) 

This can now be interpreted as the total energy of a particle moving along the direction in a 
potential U. In terms of analytical mechanics, one may declare q = ^ the canonical or Darboux 
coordinate of the motion and p = nip (p the conjugate momentum with the standard Poisson 



^ Named after Jakob Bernoulli. 
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bracket operator denoted {-,*}• The Hamiltonian of the system may then be formulated as 



J^{q,p) = l^+U{q) (56) 

Further apphcation of the chain rule of differential calculus shows that the system of equations 
obtained upon applying the time evolution operator of Hamiltonian mechanics, 

^ = -{^,5} % = -{'^^P} (57) 

is, indeed equivalent to the original equations of motion derived from the force balance using 
Newtonian mechanics. This nice relation to the Hamiltonian formalism is not apparent when 
working with the hodographic quantities u and 9. 



3.3.2 Limiting Behaviour 

Flat trajectories. A well-known approximate closed-form solution valid for quasi-horizontal 



trajectories with tan^ ^ 1 is due to Lamb (19231. The condition is equivalent to u ^ and 



necessarily requires q ^ q. Furthermore, it is required that the observed time interval be 
sufficiently short. This is clear physically, since some time or the other, the projectiles trajectory 
will turn downwards and become steeper Thus, contributions higher than first-order in t may 
be neglected. Therefore, the power series for 'P can be terminated at the quadratic term (doing 
so at the linear term itself would result in neglecting the first-order term in d#/dt), yielding 



(58) 



where ~ q for low angles of release. This yields 



1 + au.^ot 



1 + au^^ot 



(59) 



Setting V = in Eq. (14), which is tacitly assumed by both Parker (1977) and Lamb (1923) 



one obtains upon integration and elimination of t 



u. 



y,0 



+ 



(e^"^ - 1) , 



(60) 



which is the well-known solution for flat trajectories. 



Weakly resistive media. Often, one is interested in media where the effect of resistance 
is still non-negligible, but not extraordinarily pronounced, so that a may be taken to be small. In 
such cases, perturbative techniques may be used to obtain analytical approximations. One such 

2 2 

derived this rigorously by requiring tliat ^ Uy, so tliat tlie approximation may be 

made. Then the equations of motion decouple and can, in fact, be solved in closed form. Imposing the same 
condition on the solutions yields the requirements Uy o/u^ Q 1 ^^'^ ^ ^ 



Parker 



(1977 
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result is due to Tsuboi ( 1996 ), which can be rederived from the present solution. Tsuboi's solution 



is based on a first-order perturbation in the small parameter a and polynomial expansion in t; 
third-order for and fourth-order for Uy (the degree of the Uy approximate can be shown to 
be a necessary consequence of terminating the expression for at the third-order term in t). 



Therefore, the power series in the denominator of the rational expression of in Eq. (38) is 



terminated at the third-order term, neglecting quadratic and higher-order contributions in a, to 
obtain 

(61) 



This expression can now be formally manipulated using the Cauchy product to yield 



2 3 

at\-—^t — "—t + u^^QU^ 



(62) 



where again all contributions involving higher powers of a have been neglected. Finally, [Tsuboi 



(19961 made the assumption that u^ ^ 3> ^ or ttg « g (i.e. low take-off angle), so that 



,0 -atl^g'f - \Uy,,gt + {2ul,, + ul,)] . 



(63) 



Since the power series for (j) was terminated at the third-order term in the expression of u^, for 
the sake of consistency, the same should be done for Uy^ so that ^ will have to be terminated 
at the fourth-order term (recalling that = /</>dt). Using Cauchy products and neglecting 
contributions from terms higher-order in a yields 



ot 

Uy « Uy^, -gt+ t [IgH' - \Uy^,gH^ + \ {2ul^ + 3<,o) 9t - \Uy,o {'^<o + <.o)] • (64) 

Ux,0 



Expressions (63) and (64) are in agreement with the result of Tsuboi (1996) and illustrate that 



this can be seen as a special case of the general solution presented in this paper. 



3.4 Numerical Illustration 



For the following numerical examples, physical units are chosen such that a = g = 1, which 
may be achieved by chosing the terminal velocity = \fgJoL as the characteristic velocity 
and v^l g as the characteristic time scale. Also, we restrict ourselves to stationary media, i.e. 
V = or I? = w. In Fig. [T] a projectile is launched at various angles and with initial velocity 
= 1, mirroring approximately the choice of Parker (1977). For low initial angles [Q^ ^ 20°), 



the quasi- horizontal solution (60) agrees well with the numerical solution of the full problem 



obtained from a fourth-order Runge-Kutta method. For the sake of consistency, it may be noted 
that using higher-order terms in the approximation does not result in a noticeable difference. 



For a higher angle of release, such as 30°, the shortcomings of Eq. (60), which is accurate only 



up to the first order in t, as discussed in Sec. 3.3.2 become apparent. It begins to deviate from 
the actual trajectory by considerable margins. On the other hand, using higher-order terms, i.e. 



15 




Figure 1: Trajectories for the motion of a particle under quadratic drag in a stationary medium. 
The units are chosen such that a = g = 1. The initial conditions are Vg = 1 and (a) 9g = 15°, 
(b) 9o = 20°, (c) 6*0 = 30° and (d) 9^ = 40°. The 4th order Runge-Kutta solution of the initial 
value problem (16)-(19) is compared with the quasi-horizontal approximation (60) and the 
fourth-order solution obtained from the general solution (37)-(39). 



(37)-(39|, leads to a better approximation, distinctly so for 6*,, = 40°, where the assumption of a 
fairly fiat trajectory is clearly violated. 

In Fig. |2] the trajectory of an object released at Vq = 1, 9^ = 60° is depicted. Comparison is 



made between the numerical solution and the analytical solution obtained from Eqs. (38)-(39) 
upon truncation of the power series for (p at various orders n. For the sake of clarity, lower-order 
approximations are not shown in the figure. In Fig. [3| the result obtained numerically from 



solving Eq. (32) with a Runge-Kutta method is compared with the partial sums of the series 



for (p, i.e. Eq. (37). The lower order approximations are once again omitted in order to avoid 
unnecessary cluttering. From both figures, it is apparent that including higher order terms leads 
to an increase in accuracy. For a solution based on power series, this is generally a suggestive 
indication of convergence. 



4 Generalizations 

In this section two possible generalisations of the basic problem are discussed. First, more general 
drag laws are considered. In the second problem, the requirement of constant fluid velocity is 
relaxed in order to allow for the velocity profile to depend on the vertical coordinate. 

4.1 Generalised Drag Laws 

Here, the drag force is assumed to be of the form 

F^, = -nip {a\\v - Vf + p) {v - V) , (65) 
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Figure 2: Trajectory of an object released at a velocity Vq = 1 and an initial angle of 6q = 60° in 
units such that a = g = 1. The numerical solution of Eqs. (16)-(19) by a 4th order Runge-Kutta 
method (•) is compared with the solution obtained from Eqs. (38)-(39) upon truncation of ^ 
after the term of order n = 4,5, 6, 7. 




Figure 3: The resolvent variable 'P of an object released a,t Vg = 1 and initial angle 9q = 60° as a 
function of time. The numerical solution of Eqs. (16)-(19) by a 4th order Runge-Kutta method 
(•) is compared with the solution obtained from Eqs. (38)-(39) upon truncation of after the 
term of order n = 4,5, 6, 7. 
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where a, (3 and k denote real non- negative constants. The reason this form of drag presents such 
an interesting object of study is two-fold. First, from a mathematical point of view, this law is a 
formal generalisation of several simpler laws. Observe that setting a = reduces Eq. (65) to the 
well-known linear Stokesian drag. For /3 = 0, one recovers the general power law, and further 
setting k = 1/2 yields Newtonian drag. Finally, arbitrary a,/3 and k = ^ corresponds to the 
generalisation of the quadratic drag law proposed by Newton (1687). From a physical perspective. 



the form of Eq. ( 65 1 represents an often used correlation for the drag force of a sphere proposed 



by Schiller and Naumann ( 1933 ) ^[Clift et al. ( 1978 ) report that the Schiller- Naumann correlation 
is accurate to within 5% for Reynolds numbers between 0.2 to 10^. This drag law hence bridges 
the grey area between the Stokesian and the Newtonian regime. 



4.1.1 Mathematical Problem 

Inserting Eq. ( 65 ) into the general force balance ([T]) yields 

-(01117-^11"" + ^) {v-V)+g 



dt 



(67) 



Rewriting the above in terms of the relative velocity u = v — V and noting that V = const, 
yields 

[aWuf + I3)u + g. (68) 



dit 

dt 



In a Cartesian coordinate system identical to the one used in the previous sections, this equation 
may be written in components as 



du, 

At 

du, 

dt 



+ 



+ 



a{ul + uiy + p 
a {ul + ulf + a 



= 



(69) 
(70) 



The initial conditions are w(0) = {u^ Q,Uy Q) 



4.1.2 Solution 

Let us begin by recasting the equations of motion into the form considered in Section [2j setting 
P = a {ul + ul)'' + f3 and Q = (0, —g). Observe that, as in the previous problem, V„Q = 62,2 is 
still a null matrix. Therefore, this problem is of the class of equations considered in Section [2] 
and the reduction principle may be used in order to define the appropriate resolvent variable 
and derive the resolvent equation: 



exp 



a{ul + ul)'' + ^ 



dr 



^More precisely, Schiller-Naumann drag is usually stated in the form H-FuH = ICoAUd — with 



Cb(Re) = :^(l + 0.15Re°-n 
Re ^ ' 



(71) 



(66) 



where Re denotes the Reynolds number. 
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u. 



yO 



9 4>dT 



d0 
dt 



a 



"yO 



9 J^f/^d-T 



(73) 
(74) 



Introducing <P according to Eq. (31) yields an ordinary differential equation of second order: 



a 



ulo + {Uyo - 9^) 



fe /d^ 
d^ 



d^ 

d^ 



(75) 



The initial conditions remain the same as in the previous problem. A power series ansatz is now 
used, setting 



^ = 6^ + 6^^ + h^e + h^r" + 



(76) 



The initial conditions require that 6o = ^I'^d h^ = 1. All further coefficients then follow from 
Taylor's theorem and may be computed recursively using 



a d" 



a 



(n + 2)! dr 
The first few terms of the solution hence read 

6n = 



ulo + {Uyo - 9^) 



~dt 



+ /5 



,d^> 
d^ 



(77) 



6, = 1 



(1 - 2k) a'uf - 2k aguf-^ sin + 2 (1 - A:) a^uf + 

2^~3 



Inserting this solution into Eqs. (72) and (73) completes the solution of the problem. 



Remark. The general nature of the coefficients is a bit more transparent if one chooses 
units such that = g = 1. Then the coefficients are bivariate polynomials in a and (3. Each 
coefficient comprises terms of three different kinds: First, there are those terms that depend only 
on a and thus represent the effect of the general power drag law, which is associated with high 
Reynolds numbers. Secondly, there are terms that represent the effect of Stokes drag, i.e. low 
Reynolds numbers. Other than these pure terms, there are mixed terms that are responsible 
for the transition between these two regimes. Although the underlying drag law is, formally, a 
linear superposition of a general power law and a linear law, the solution itself is not a sum of 
the solutions corresponding to the "pure" drag laws. As such, the existence of the mixed terms 
illustrates the nonlinearity of the problem. Nonetheless, this behaviour only becomes apparent 
from the third-order term onwards. 
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It is readily checked that upon setting /3 = and k = 1/2, one obtains the same expansion as 
estabHshed in Section 2. On the other hand, for a = 0, one obtains 'P = t+ f3t^ /2\ + f3^t^ /3l + ■ ■ ■ = 
(e''* — 1) //3, whence = u^^/e^^^ and Uy = (n^ g — g/l3)e~^^ + g/ (i-, which is the weh-known 
closed-form solution for Stokes drag (cf. Synge & Griffith 1949, p. 159). This serves as a 
consistency check. 

4.2 Accounting for a Velocity Profile 

As a final step of generalisation, we consider the effect of velocity profiles on projectile motion. 
After deriving the mathematical formulation, i.e. the equations of motion, the solution is obtained 
using the result of Section |2j Finally, some properties of the solution are pointed out. 

4.2.1 Equations of Motion 

The equations of motion for projectile motion subject to quadratic drag are given by 

— = -a\\v-V\\{v-V)+g. (78) 

The usual Cartesian coordinate system is chosen. Unlike in the previous problems, however, the 
fluid velocity is now allowed to vary with height. We assume a unidirectional horizontal flow 
along the x axis with the velocity depending only on the vertical coordinate y, i.e. 

^ = (K(y),o) (79) 

so that 14 is a sufficiently well-behaved function of y, subject to the condition V^{y = 0) = 0. 
Physically, this represents a flow in the x-direction over a horizontal plate or wall placed at 
y = so that the fluid obeys the no-slip condition of a viscous fluid at a solid wall. Here, we 
restrict ourselves to the case where the velocity of the fluid varies linearly with height, setting 
= lU- Such a law is valid for flows (both laminar and turbulent) close to a smooth wall. Now 
the differential equations describing the motion may be written in coordinate form as 

'^ + a{v^- 7y) ^(v^--iyf + vl = (80) 
+ avy^ {v^ - lyf + vl = -g (81) 

The initial conditions are, as usual, v^{0) = v^ g and Vy{0) = Vy ^. In addition, it is assumed that 
the particle is released into the flow at x{d) = and y(0) = y^. Recalling that Vy = dy/dt, one 
recognizes that the equations of motion actually constitute, in their present form, a second-order 
system of coupled differential equations (the second equation involving dVy/dt = d^y/di^) and 
are, as such, beyond the scope of the reduction principle introduced earlier. This inconvenience 
can be resolved if the problem is reformulated in terms of the relative velocity u = v — V = 



20 



- ly, Vy) = {u^, Uy) yielding 



dt 
du 

~dt 



" + au„\/ul + ul = -g 



(82) 
(83) 



subject to the initial conditions n^(0) = u^ q and Uyifi) 



u. 



y,o- 



4.2.2 Solution 



In order to recast the system into the form considered in Eq. the formulation of P may be 



adopted without change from Eq. (26), since the same form of drag applies here as well. On the 



other hand, as an effect of the velocity profile, the inhomogenous term now reads Q = {—jVy, —g). 
The gradient of Q hence is 

.0 ' 



(84) 



i.e. no longer zero. Therefore, the traditional methods such as the ones used by Bernoulli (1719) 
or 



Parker (1977) are inadequate. However, the problem is still within the scope of the reduction 



principle from Section [2j because V^Q is, up to numeration of indices, a strictly lower triangular 
matrix. The procedure then yields 



u,. 



exp <; a / \l ul + ul dr 



Uyo- 9 [ <pd.T 
Jo 



UxQ-1 Uy{(l>)<Pdt 

Jo 



w.o - 7 / [Uyo-9 0ds dr 



(85) 
(86) 
(87) 



This may be further simplified by setting 



0/ = 



)dMt. 



The first initial condition is imposed by the definition of i;^ = d^'I^/dt^, namely 

d^<f (0) _ ^ 



df 



(90) 



Two other initial values are still necessary in order to determine !P' uniquely. This may be 



accomplished, e.g., by fixing the choice of integration constants in Eq. (89). Keeping in mind 



21 



^(0) = 




91 


d<P'(0) 


—u 


dt 


9 



that the physical initial conditions u^ifi) = u^ ^ and My(0) = Uy ^ need to be fulfilled, we choose 

(91) 
(92) 

Now, the solutions may be rewritten in terms of the auxiliary function ^ as 

dip-/dt , , 



Inserting this into Eq. ( 83 ) , yields 



At' 



agJ{j^) + {d>^/dty. (95) 



This is the new auxiliary equation to which the original coupled system has been reduced. One 
can now develop ^P' in a power series = Cq + c-^^t + Cat^ + c^t' + • • • with 

c„ = -r-^ (96 
" nl dt" ^ ' 

according to Taylor's theorem. The first three coefficients are given by the initial conditions. 



Differentiation of Eq. ( 95 ) n times yields the recursive formula 

which can be used to evaluate all the coefficients of the series solution, the first few terms reading 

Uq cos ^0 

Co = 

91 

Un sin On 



9 

c, = 1/2 

^0 



aun 



2-3 

a {g sin 6*0 — iUq cos 6^ sin Bq) 
2 • 3 -4 
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The solution then is 



ag sm (/q — aj sm cos „ 
1 + auot 1 + 



1,2 5a^^o ,3 , a (5' sin 6*0 -57 cos sin 00 ) 4 (98) 

UnSm&nt 1 1 H 1 +••• 

X ° ° 2 1-2-3 1 •2-3-4 

ag( sin 6*0 — 07 sin cos „ 
1 + augt 1 H 



ag( sm u„ — 07 sm cos „ 
1 + augt ^— ^ 1 + 



1 + + g (g' sin 6*0 - 57 cos sin Oq) (99) 

X 1-2 1-2-3 

ag sin 6*0 — ju^ cos sin 
1 + auot — \ 

Using the transformation equations given earher, one can deduce v^,Vy,x and y from the above. 



4.2.3 Properties 



Comment on the auxiliary variable. Here, in Section 4.2 the abstractly defined resolvent 



variable no longer has the usual connotation of being (up to linear transformations) the slope 



of the trajectory of the projectile. This also explains why the usual approach (cf. Parker 1977) 
of assuming a resolvent variable linearly related to Uy/u^ is not fruitful in this case. It may be of 
interest to note that 

= (100) 

Ux 7 5^ 

which is, up to constant factors, the logarithmic derivative of iP". From a similar point of view, 
the hodograph technique also appears to be problematic because expressing 'F in terms of 



tan^ = Uy/u^ and inserting the resultant expression into Eq. (95) would yield an integro- 
differential equation involving non-trivial kernels. 

Particular Solutions. Although the complete analytical solution could not be expressed in 
closed form, surprisingly, there are several distinct non-trivial particular solutions that may be 
expressed in terms of elementary functions. To this end, assume an exponential ansatz of the 



form ^ = Cexp(At), where C is an arbitrary constant. Inserting this into Eq. (95) yields the 
characteristic equation 



A3 = 0(^^72 ^A^ (101) 

for A. This equation has three distinct roots in C corresponding to three separate solutions for 
the set of initial conditions ^{Q) = C,dl^(0)/dt = AC, d^!P'(0)/dt^ = A^C. Since there is only one 
free parameter contained in the solution, a general set of initial conditions cannot be fulfilled, 
neither can the general solution be obtained by linear superposition owing to the nonlinearity of 
the problem. It would appear, nevertheless, that the properties of the particular solutions reflect 
the behaviour of the general solution, so that knowledge of the former would be instructive in 
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2 4 6 8 10 



Figure 4: Family of phase trajectories for in normalized units such that a = g = 1 and fixed 
7 = 1. The initial conditions are varied with u^ g = 1 and g = 0.5, 1, 1.5, and 2. 



comprehending the latter. This may be visualized by means of phase plots in the {^,d^/dt) 
plane. In Fig.|4] a family of phase trajectories is displayed for 7 = 1 and various initial conditions, 
with a = g = 1 achieved by an appropriate choice of units. The qualitative shape may indeed be 
derived from the nature of the particular solutions named before. Using elementary algebra, it is 



readily seen that the solutions of Eq. (101 ) are such that one is positive real and two of them 
are complex conjugates with negative real part. The complex conjugate solutions are associated 
with exponentially damped sinusoidals, translating into a spiralling ellipse in the phase plane. 
Their influence, however, decays with time due to the damping. The real root, representing a 
growing exponential function, quickly begins to dominate the solution, so that the long-term 
behaviour of the solution in Fig. |4]is a straight line. The slope is given by the real solution of 
= VFTl, i.e. 



A ' 

^3 




1.15096 



5 Conclusions and Outlook 

In this paper, the motion of a body in a resistant medium and a constant gravitational field 
was investigated. The nonlinear nature of the drag force leads to a coupled nonlinear system 
of ordinary differential equations, thus demanding a more involved approach in the analytical 
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Figure 5: Situation of a particle released at y = 0.05 in units such that a = g = 1 and with 
initial conditions v^ g = l,Vy Q = 1/5. Displayed are the trajectories for different values of the 
slope of the velocity profile, 7. The analytical solution obtained from truncating \^ after the 
fifth-order term (solid line) is compared with the numerical solution by a Runge-Kutta method. 
For reference, the Runge-Kutta solution without velocity profile (7 = 0) is also depicted (dotted 
line) . 
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treatment. To this end, a useful technique for decouphng a particular class of coupled systems of 
ODEs was proposed in order to facilitate the analytical solution of such problems. 

This was applied to projectile motion under quadratic drag, yielding an equation of motion 
for a single auxiliary variable. An explicit analytical solution of the problem was proposed. 
Relations to limiting cases from earlier studies were then shown. The solution was also validated 
numerically by comparing it with direct solutions of the original system of equations using 
a Runge-Kutta method. Furthermore, although a constant of motion for this problem was 
already known in the literature, it was further demonstrated that, when expressed in terms of 
appropriate variables, the constant of motion can be interpreted as the conserved Hamiltonian 
of the system. Within this framework, it was further pointed out that the resolvent variable 
provided by the reduction principle together with its time derivative constitute the canonical or 
Darboux coordinates of the motion. 

The decoupling technique was carried over to handle a more general drag law containing 
quadratic and general power drag laws. A special case of this law is the Schiller-Naumann 
correlation valid for the transition regime between Stokesian and Newtonian drag. An explicit 
analytical solution was presented for this well. 

Lastly, the reduction principle was also used in order to deal with a final generalisation 
of the basic problem, the motion of an object with quadratic drag in shear flow with a linear 
velocity proflle. This was motivated by the fact that no study investigating projectile motion 
under the effect of a velocity profile and nonlinear drag could be found in the literature. Again, 
an analytical solution was presented. Furthermore, closed- form nontrivial particular solutions 
for the resolvent variable were identified and the close relation between their properties and the 
nature of the general solution was discussed. These could be the basis of new studies, either 
from the perspective of nonlinear superposition principles (i.e. with the aim of constructing 
a general solution from the particular solutions) or from the point of view of semi-analytical 
approximations, e.g., by means of the method proposed by Churchill and Usagi (1972 1 . 

While in practical applications a numerical solution of the class of systems of ODEs considered 
here may be more efficient in computing a solution for given initial conditions, the present 
results are important as they provide a means of access to the mathematical properties of and, 
ultimately, more physical insight into the system, even in more general cases, as demonstrated. 

Several extensions of the present work are now possible: 

The proposed decoupling technique may be applied in other settings. These need not be 
within the direct scope of physics; coupled systems also occur in the modeling of chemical or 
biological reactions, to name but few. An extension of the method to other possibly more general 
classes of differential equations may also be of interest. 

Second, it seems appealing to generalise the shape of the background velocity profile of the 
fiuid to some classes of nonlinear functions. Finally, the considered problem can be enhanced by 
considering other forces as well, of which the Magnus force presumably is the most significant, 
requiring an additional equation for the angular momentum of the particle. 
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